f(R) cosmology revisited 



Luisa G. Jaime 1, 2 Q Leonardo Patino 2 Q and Marcelo Salgado^ 
1 Instituto de Ciencias Nucleares, Universidad National Autonoma de Mexico, A. P. 70-543, Mexico D.F. 04-510, Mexico 
2 Facultad de Ciencias, Universidad National Autonoma de Mexico, A. P. 50-542, Mexico D.F. 04510, Mexico 

(Dated: March 7, 2013) 

We consider a class of metric f(R) modified gravity theories, analyze them in the context of 
a Fricdmann-Robertson-Walker cosmology and confront the results with some of the known con- 
straints imposed by observations. In particular, we focus in correctly reproducing the matter and 
effective cosmological constant eras, the age of the Universe, and supernovae data. Our analysis 
differs in many respects from previous studies. First, we avoid any transformation to a scalar-tensor 
theory in order to be exempted of any potential pathologies (e.g. multivalued scalar potentials) 
and also to evade any unnecessary discussion regarding frames (i.e. Einstein .vs. Jordan). Second, 
based on a robust approach, we recast the cosmology equations as an initial value problem subject 
to a modified Hamiltonian constraint. Third, we solve the equations numerically where the Ricci 
scalar itself is one of the variables, and use the constraint equation to monitor the accuracy of the 
solutions. We compute the "equation of state" (EOS) associated with the modifications of gravity 
using several inequivalent definitions that have been proposed in the past and analyze it in detail. 
We argue that one of these definitions has the best features. In particular, we present the EOS 
around the so called "phantom divide" boundary and compare it with previous findings. 
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I. INTRODUCTION 



Astronomical observations based on type la supernovae (SNIa) together with the assumption that the Universe 
is homogeneous and isotropic at large scales led to the conclusion that the Universe is currently expanding in an 
accelerated way [3-[3] ■ This phenomenon can be most easily explained by appealing to the existence of a cosmological 
constant A (sometimes termed dark energy) . This constant along with the introduction of dark matter (DM) appar- 
ently needed in many regions of the Universe (galaxies and clusters) have originated what is called today the ACDM 
paradigm. This paradigm has also successfully explained most of the current details of the Cosmic Background Ra- 
diation (CBR or CMB) in the framework of general relativity as well as other important features of the Universe 
at large scales [5| (for a thorough review see Ref. |6|). 

However, despite of the simplicity and success of this paradigm, several theoretical as well as epistemological 
arguments have been put forward as objections against such a simple model of the Universe. For instance, as concerns 
the DM hypothesis, one of the the main criticisms is that its nature (i.e. its quantum and classical properties) is not well 
understood (if at all) yet. That is, apart from the gravitational evidence, there is no further strong reason supporting 
its existence. Since several experiments have failed so far to detect the proposed DM particles, skepticism keeps 
growing in this direction. On the other hand, the cosmological constant has been historically regarded as "suspicious" 
by several detractors (including Einstein himself; see Refs. Q for a review), although some of its apparent drawbacks 
are based more on prejudices than on strong and well grounded physical arguments (sj. In any case, the discomfort 
that A has produced in the spirit of some people has led to consider more complicated alternatives, of which, is fair 
to say, none is regarded today as a more serious candidate for dark energy (DE) than A (the BigBOSS experiment Q 
has been designed to shed light in this direction). Among these alternatives are the so called modified theories of 
gravity (MTG) as opposed to general relativity (GR). Some of these theories have been also proposed to substitute 
DM and even as models for inflation. Perhaps the most popular MTG over the past ten years and the one we focus 
in this article are f(R) metric theories, where an a priori arbitrary function of the Ricci scalar R replaces R itself in 
the gravitational Lagrangian. 

This kind of MTG were conceived originally in order to create a late accelerated effect without a cosmological con- 
stant or as inflationary model without an extra scalar field (see Refs. |10l4l3l| for a detailed review). Notwithstanding, 
despite of the promising f(R) models first proposed to replace A 14], a cumulative evidence, both theoretical and 
observational, has been found against most of them. However, new models have been proposed to overcome the initial 
difficulties, some better motivated that others but none introducing a new fundamental principle that can be used as a 
guiding line; indeed they have rather been constructed by trial and error. The simplest (non-trivial) choice f{R) = R 
was historically favored by Einstein since mathematically led to second order partial differential equations (PDE's) 
which could easily reduce to the Newtonian theory in the week field limit. 

General mathematical and physical conditions are usually demanded in order to avoid pathologies in the models. 
For instance, the conditions /as > and fa > (where the subindex indicate derivative with respect to R) seem to 
be required for stability considerations and to ensure a positive definite effective gravitational constant, respectively. 
It is however not clear if those conditions are really necessary, and in many studies they are not imposed. Therefore, 
in most cases, "handcraft" has been used to design a particular f(R) model based on heuristic arguments that might 
account for the actual phenomenology when the full fledge model is submitted to a detailed scrutiny. The general 
trend so far is that no single model is able to explain most of the current observations, but only some aspects of them. 
That is, most of the f(R) models fail miserably when they are preempted as models for all the dark substance (both 
DE and DM) and when taking into account the Solar System tests as well. Even when considered only as DE models, 
most of them fail, with the exception of some notable cases. Of course it could well happen (but perhaps not very 
desirable) that dark matter, dark energy and modifications of the laws of gravity in some combination are required 
by nature in order to fully understand our Universe. In fact, this is the approach we pursue here at the cosmological 
level except that we do not include explicitly a cosmological constant, but rather, the models themselves give rise to 
an effective A which vary slightly around the required value at late (i.e. present) epochs of the Universe. We then 
consider f(R) theories simply as a model for explaining the accelerated expansion of the Universe and introduce a 
dark matter component in the same way as in the ACDM paradigm. Nevertheless, things turn out to be not so 
simple, given that the proposed models can disturb the successes of GR. Any proposed specific f(R) model has not 
only to satisfy the cosmological observations, but all the gravitational observations at all scales. These include the 
Solar System experiments, the existence of physically acceptable compact objects, the binary pulsar, etc. As today, 
there is no single MTG model that replaces successfully GR and explains as should be, all the observations for which 
it was designed originally. 

After the first f(R) models were proposed to explain the accelerated expansion of the Universe, among them the 
"historic" f(R) = R — /i /R, a sequence of papers appeared where the constraints imposed by the Solar System were 
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taken into account [15j . Without reaching a clear consensus on the issue, it seemed that such models were not viable. 
One of the arguments put forward to establish that conclusion was based on the fact that such theories can be shown 
to be dynamically equivalent to a Brans-Dicke (BD) theory with ui = 0. Since such a value for uj in these theories 
gives rise to a post-Newtonian parameter 7 = 1/2, which conflict with 7 ~ 1 favored by the Solar System tests, then 
at first sight the analysis suggested that all f(R) theories were excluded blatantly as viable theories. Much later, 
it was recognized that such an argument should be used with care in view that f(R) theories are not equivalent to 
the standard BD theory with to = 0, but to a BD theory with a potential. Therefore, depending on the mass of 
the effective scalar, the theory at hand could pass or fail the Solar System tests (lq . Although it is now recognized 
that many of the f(R) models give rise to 7 ~ 1/2, and are therefore ruled out, some others, due to the effective 
mass of the scalar, might produce a successful phenomenology. This success depends on whether or not the scalar 
field which is associated with the model at hand can act as a chameleon [l7l . fl8j . a mechanism that appears in 
some scalar-tensor theories of gravity which allows them to satisfy the local tests and the possibility of producing the 
required cosmological effects . 

Now, as concerns the cosmological restrictions on these theories, Amendola et al. 0, [2l[ have devised criteria of 
quite general applicability that allows to discard many of the proposed f(R) models. In short, their analysis shows 
for a large class of models that they either produce an accelerated expansion at recent times but fail to generate a 
correct matter-dominated era (the scale factor behaves as radiation in GR) or the opposite. In fact, when viewed 
from the past to the present, many of such models cross from the radiation-dominated era (deceleration epoch) to 
an effective dark-energy era (acceleration epoch) with a very short or not even existent matter-domination era. Such 
models are therefore incompatible with the CBR observations, the age of the Universe and the structure formation 
at large scales [13, [H|. This issue was, however, not free of debate either [22I l23l|. 

To make things even more confusing in this matter, additional skepticism was raised about the viability of such 
kind of theories when a further test was performed on several cosmologically successful f(R) models. This time, the 
test consisted in analyzing the possibility that such models allowed the construction of solutions representing realistic 
(or at least idealized) neutron stars. In a first attempt to do so, Kobayashi & Maeda [24| showed that idealized 
neutron stars (specifically, incompressible compact objects) were not allowed by the Starobinsky model [25j since a 
singularity in the Ricci scalar developed within the object. Later, Babichev & Langlois 26| criticized such conclusion 
and argued that the singularity was only due to the use of an incompressible fluid and when a similar situation was 
analyzed with a compressible gas (e.g. a poly trope) no singularity was found. Finally, Upadhye & Hu [13] argued 
that a chameleon effect was the responsible of avoiding the formation of singularities in compact objects and not the 
use of more realistic equation of state. 

In a more recent work by us [28| , while we arrived to the same primary conclusion of Refs. (26l [27j on that no 
singularities were necessarily formed, we criticized some very basic aspects of the analysis common to all of the three 
above investigations concerning the existence of neutron stars [2J, |26| . |27| . To be specific here let us mention that 
their analysis relied fundamentally on the fact that a potential associated with a scalar-tensor counterpart of the 
Starobinsky model could drive the scalar field to a point where the corresponding Ricci scalar diverged. The main 
point of our criticism, noticed earlier in [29l [30|. is that such a potential has unpleasant features as it is multivalued 
and therefore, the transformation to the STT is not well defined. Since the conclusions reached on those papers 
depend crucially on the use of such a pathological potential we consider them not very trustworthy even if the 
dynamics is supposed to take place in a region where the potential is single valued. Moreover, due to the fact that 
many of the controversies regarding this subject have been the result of using the STT in the Einstein or Jordan 
frames (see also Ref. [3l| for a similar criticism), we strongly suggested to abandon such an approach and treat 
f(R) models in all applications without performing the transformation to any frame of STT. In the particular case 
of static and spherically symmetric spacetimes, we devised a very transparent and simple approach that allowed us 
to deal with compact objects. Using a particular case of the Starobinsky model we concluded that no singularities 
were found. We also studied the Miranda et al. model [3(| (which was previously shown to be free of singularities 
following the STT approach but with a single- valued potential), and arrived to the same conclusion. We emphasized 
the advantages of using that robust approach over the STT technique, and we argued that even in the cases where the 
STT transformation is well defined it is not particularly useful and that the original variables are required anyway in 
order to interpret the solutions correctly. Furthermore, and as we mentioned above, dealing with the STT approach 
opens the way to the long-standing controversy about which frame (Einstein or Jordan) is the physical one [321 ] . This 
discussion is in some cases semantic and in many others completely ill founded and corrupted. In our treatment we 
don't even need to deal with it at all since we consider the theory directly as it emerges from the original action without 
performing any formal or "rigorous" identification or transformation with any other theory, frame or variables. 

In the present article we review, in light of our robust approach, the cosmological analysis of some of the apparently 
least problematic f(R) models, in the sense that they seem to pass the cosmological and Solar System tests (although 
the model of Ref. [30( is still on debate). As we will show, our treatment allows to handle the equations most as 
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in the case of GR, which in turn, makes possible the use of simple techniques of numerical relativity to monitor the 
accuracy of the solutions. 

The paper is organized as follows, Section II introduces the f(R) theory, the general equations and the approach 
under which we will treat them. In Section III we present four inequivalent definitions of the "energy-momentum" 
tensor of geometric dark energy that we have identified in the literature and which give rise to three inequivalent 
equations of state (EOS) used in cosmology. Section IV displays in detail the three specific f(R) models that we 
submit to a cosmological analysis. In Sections V and VI we focus on the Friedmann-Robertson- Walker (FRW) 
cosmology and analyze the EOS that arises in the particular f(R) models we treat, as well as the relative abundances 
of the different matter-energy components in order to explicitly show the matter and modified-gravity (geometric dark 
energy) domination epochs. In order to have some insight about the viability of the solutions, we compare the results 
with the successful GK-ACDM scenario. We also compute the luminosity distance and confront the results with the 
historic SNIa data Q and the UNION 2 compilation Q. The age of the Universe that arises from these models is also 
estimated. Finally, Section VII concludes with a summary and a discussion. An Appendix displays the dimensionless 
form of the cosmological equations used for numerical integration and the numerical test used to check the accuracy 
of our solutions. 

II. f(R) THEORIES, A ROBUST APPROACH 

The MTG that we consider is given by the following action 

S[g ab ,il>] =J^V=9 d 4 x + S niatt {g ab ,i>] , (1) 

where K = 8ttGo (we use units where c = 1), and f(R) is an apriori arbitrary function of the Ricci scalar The 
first term corresponds to the modified gravity action, while the second is the usual action for the matter, where ip 
represents schematically the matter fields (including both the visible and possibly the dark matter). 
The field equation arising from Eq. (fTJ) in the metric approach is 

fnRab ~ T^fgab ~ (V a V& - #a6 D ) /fl = K,T ab , (2) 

where f R indicates d R f , □ = <? afl V a Vb is the covariant D'Alambertian and T a b is the energy-momentum tensor of 
matter which arises from the variation of the matter action in Eq. (TTJ). It is straightforward to write the above equation 
in the following way 



~ (RfR -f)+ fRR^R + f R R R (VRj 2 



KTab , (3) 



f R G ah ~ f R RV a V b R - fRRR ( V a R) ( V&iJ) + g ab 

where (Vi?) 2 := g ab (V a i?)(Vfci?). Taking the trace of this equation yields 

kT - 3f RRR (VR) 2 + 2/ - Rf R , (4) 



UR= 1 



3Jrr 

where T := T a a . Finally, using Eq. (gj) in Eq. © we find 



G a b — -7— 

JR. 



/flflV a V b i? + f R R R {V a R)(V b R) - ^ (Rf R + f + 2 k t) + KT ab \ . (5) 



Equations Q and §5§ are the basic equations for f(R) theories of gravity that we propose, as previously stated 
in (28|, to treat in every application, instead of transforming them to STT. Notice that GR with A is recovered for 
f(R) = R — 2 A. It is important to stress that even in the treatments where the STT approach is not pursued, in 
most of the cases the term DR is not rewritten in the way we do it here, and therefore the resulting equations for 



1 It is important not to confuse this parametrization in the action with the alternative writing f(R) = R + f(R) which is used by several 
authors (c.f. [iHl), and where the tilde is then dropped. 
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specific spacetimes turn out to be much more involved. The idea is also to rewrite, when possible, all the second order 
derivatives of R coming from V a Vfci? in Eq. ([S]) in terms of lower order derivatives using Eq. We shall do that 
for the cosmological applications that we analyze in Sec. [V] The systematic approach to perform the full first order 



reduction is through the 3+1 formalism 331 ] . The method proposed here follows the same philosophy of our previous 
article [28| . We also point out that a similar reduction was considered by Seifert t 34] , when analyzing the stability of 
several systems under the framework of f(R) theories. 

An important property of these theories is the well known fact that not only the total "energy-momentum" tensor 
(i.e. the right-hand-side -r.h.s- term of Eq. [5]) is conserved, but also the energy-momentum tensor of matter alone 
T ab . That is, the field equations imply V a T ab — 0. This reflects no other but the fact that this kind of theories are 
metric theories, and so, the geodesic equation for test particles holds. The matter equations will take then the same 
form as in GR, and the departure from the latter will occur in the evolution equations that the metric will follow. 

Clearly, we have not modified the fourth-order character of the theory with respect to the metric, since R depends 
on second derivatives of the g a b components and, as we appreciate from the r.h.s of Eq. (0, there are in addition 
second derivatives acting on R. Nevertheless, the important point here is to promote R as an independent variable, 
and solve the system @ and ([5]) as a set of coupled second-order PDE's for R and g a b respectively. Of course, the 
same spirit is used when the theory is transformed to an STT counterpart, where instead of R a scalar-field x = fn 
is defined. However, as we emphasized before, we will avoid such a treatment in view of the potential drawbacks that 
can appear, while in our case everything is as well defined as the function /(i?) itself. In particular, since we shall 
deal with the Starobinsky [25[ and Hu-Sawicky [l8[ (hereafter HS1) models where fun is not positive definite and 
for which the scalar-field potential that arises in the STT transformation is multivalued, it is then advisable not to 
pursue that approach. 

It is important to stress that following this treatment (i.e. the second-order approach as opposed to the fourth 
order one for the metric alone), the initial data on R and its time derivative are not arbitrary, but subject also to a 
modified Hamiltonian and momentum constraints [33]. These important mathematical issues become apparent when 
formulating the theory as an initial value problem [33|]. In the particular case of a FRW spacetime, this issue will be 
reflected in that the initial values for R and R must satisfy the equivalent of the Friedman equation for GR, which 
amounts to the modified Hamiltonian constraint; the momentum constraint being trivially satisfied in this case. Under 
this approach, there exists another identity which links the Ricci scalar with first and second order derivatives of the 
metric. This identity can provide a consistency test for the numerical integration since, like the initial data constraints, 
it must be satisfied everywhere in the space-time. In particular, in static situations where there is no evolution or 
the evolution is trivial, this identity can be used to check the self-consistency during the numerical integration over 
the space [28]. A similar situation happens for the FRW case, except that in this case R provides another evolution 
equation for the Hubble expansion. Nevertheless, this equation is consistent with the rest, and so, it provides in 
practice, a redundant equation, since a priori no extra information can be extracted from it, although it can be used 
also to fix the initial data (by means of the deceleration parameter and the jerk - see Sec. IVBp . As in the static and 
spherically symmetric case, this redundancy can be exploited numerically to check the consistency and accuracy of 
the computer codes that we created to solve the equations numerically. 

Finally, we mention what maybe the most important property of f(R) theories as models intending to mimic a 
cosmological constant at present time. First, from Eq. (jj) a "potential" V(R) = —Rf(R)/3+ J R f(x)dx can be 
defined, so that Vr(R) = dV(R)/dR = (2/ — fuR) /3. Then notice that if the matter terms are absent (e.g. outside 
a compact object), almost negligible or small (e.g. at late times of the Universe), then Eq. @ admits as solution 
R = Ri = const., provided Vr(Ri) = (c.f. Sec. II VI for specific examples). When this solution is used in (JSJ) under the 
approximation T a b ~ 0, this reads G a b = —-^eS9abi where A e ff := i?i/4. Of course taking the trace of (J5j) is consistent 
with R = R\. Then the solution for the metric must be such that the Ricci scalar is constant. But this is no other 
than the de Sitter type of solutions. In particular, this holds for static and spherically symmetric spacetimes [HI [H[ 
and for the FRW cosmology, as we shall see below. So, in summary, f(R) are able to mimic A provided that one 
finds solutions where the matter contribution can be neglected, e.g. asymptotically in space or in time, and where 
the Ricci scalar R approaches a critical point (a maximum or minimum) of the "potential" V(R). All this relies on 
Jrr(Ri) 0, as otherwise the conclusions might change. Given that we solve the full equations ^ and (O regardless 
of any ad hoc definition of a potential, in our case V(R) is merely used as a guiding tool for identifying the critical 
point associated with the possible asymptotic solution for R. Then, only for such guiding purposes and for models 
where fan is positive definite in general or that < fan < oo in the regions where the solutions take place, it is 
irrelevant to include the term fan in V(R), since we will find the same critical points. Nevertheless, when those 
conditions are not fulfilled (i.e. when Jrr{R\) — or when Jrr{R\) — > oo), it is then advisable to include the term 
fun in the potential, otherwise one could "miss" a critical point (see Sec. HVI for a further discussion). For the f(R) 
models that we analyze below, some of which have a non positive definite fuR (see Sec. IIV[) it turns however, that 
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the cosmological solutions presented in Sec. IVI1 never reach the point(s) where Jrr < nor where Jrr diverges. 
Therefore, there exists viable cosmological models that are free of any "pathologies" of this sort. 



III. THE "ENERGY-MOMENTUM TENSOR" OF f(R) 



Very often it turns to be convenient and helpful to write the field equations of alternative (metric) theories of 
gravity as the Einstein field equations with an effective (total) energy-momentum tensor (EMT) that contains all the 
modifications which are associated with the new theory and the EMT of matter itself. In cosmology this rearrangement 
of the equations can be useful as one tries to identify the contributions of the modifications of gravity within the total 
EMT as though they represented some kind of geometric dark energy, and so it will be dubbed. This can be specially 
advantageous since one can define then an EOS associated with such dark energy and compare it with the K.CDM 
model. The only problem with this construction is that, given such total EMT, there is no canonical way to perform 
the separation between the matter and the geometric dark energy contribution, and thus, very often different authors 
introduce different definitions for the EMT of geometric dark energy. The reason is that even the matter energy- 
momentum tensor in Eq. ([5]) appears to be multiplied by a factor fZ , and so, it would seem a priori difficult to 
unambiguously define which part of the field equations belongs purely to the matter terms and which corresponds to 
the geometry (c.f. the remarks at the end of Sec. IIA of Ref. (12l|). 

In the following we provide four inequivalent definitions of the EMT of geometric dark energy, that when applied 
to cosmology, lead to three inequivalent ways of defining its corresponding EOS. In our opinion, the fact that these 
various definitions have been considered in the literature without even identifying them as inequivalent, has added a 
great amount of confusion to the subject. We hope that by clearly exposing these definitions we help to clarify this 
matter. 

Recipe I: First, define kT" as the r.h.s. of Eq. ([5]). Then, define the EMT of the geometric dark energy as 



rpX 

^ab 



1 ab 



Tab 0- By construction, T^l is conserved, for T*?* is conserved by the Bianchi identities, and 



as 



mentioned in Sec. II, the energy-momentum tensor T ab of matter alone turns to be also conserved (c.f. Appendix A 



1 ab 



T a b, which in turn leads to T* b 



of Ref. |35| for similar considerations and also Ref. [25| for further reflexions). Moreover, in the GR case f(R) = R, 



0, even in the presence of matter, as one can verify from Eq. © below. We 
conclude that includes a non trivial contribution only when GR is modified, and thus captures the idea about the 
energy-momentum content of the geometric dark energy. The explicit form of T^ is as follows: 



T 



x 



1 



fRR^aVbR + f R R R {V a R){V R) - ^ (RfR + / + 2 K T) + KT ab (I - f R ) 



(6) 



As we mentioned above, clearly for f(R) = R, T^ = 0, and for f(R) = R — 2A, — —kg ab /n, where we have used 
R = -nT tot . 

It may perhaps seem awkward to see the EMT of matter T ab appearing in the definition of the EMT of geometric 
dark energy Eq. ([B]). Nonetheless there is a simple way to rewrite T^ h in terms of purely geometric quantities. Adding 
G a b to both sides of Eq. ^ and defining T* h such that G ab = kT^ + nT ab , we obtain 



KT^ = G ab (l - f R ) + /flflVaVfci? + fRRR(V a R)(V b R) - 3 ab 



Wr -f) + fRR°R + fRRR.(VRf 



(7) 



Alternatively we can replace OR using Eq. and get 



kT* = G ab {l - f R ) + f RR V a V b R - 



f RRR (V a R)(V b R) - 9 -f (rJr + f + 2kT 



(8) 



the 



which of course can be recovered by using nT a b — G a b — K T^ in Eq. ([6]) and then solving for . Like in Eq. 
matter contribution appears also in Eq. (JSJ via the trace T = T a a of the EMT of matter. 

We have then arrived to three equivalent ways of writing T^j which are given by Eqs. ([6]) — ([8]), except that now 
in Eq. ([7]) the matter contribution does not appear explicitly. Expression ([7]) corresponds exactly to the EMT of 
geometric dark energy considered in Refs. [25|, |3&|38| 0. Following the line of thought that we have proposed in Sec. II, 



2 In several articles is denoted by instead. Nevertheless, sometimes it goes beyond purely notation and different meanings are to 
be understood in both symbols (c.f. Recipe III). 

3 In Ref. [25l| the signature (+,—,—,—) is used along with a different sign convention for its TR E . 
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and for the cosmological applications, it will be better for us to work with Eq. ([6]) rather than with the other two 
equivalent expressions. 



Recipe II: In order to obtain the second proposal for the EMT we follow a prescription similar to Recipe I, but 
we write the EMT of the geometric dark energy as the following linear combination T* b ' X (A) := AT^ — T ab — 
An~ l G ab — T ab where A is a constant. This is a rather ad-hoc generalization of which in our opinion has no deep 
motivation. By the same arguments given above, this EMT will also be conserved. Explicitly it reads 



kT%< x (A) :-- 



G ab (A - f R ) + f RR V a V b R + f RRR (V a R)(V b R) - g ab 



Wr -f) + fRR°R + fRRR^R? 



, (9) 



so that Eq. © reads AG ab = kT'JJ* + nT ab . Thus, = T ab ( 1 ) , that is, T^ b ' and coincide if and only 
if A = 1. Several authors |2ll . I39l44l| considered this recipe in the cosmological context and set A = f R (Fq in their 
notation; the knot indicating today's value). Therefore when f R is not taken as unit, as may be often the case, the 
EOS obtained from Recipes I and II are not equivalent (see Sees. IV M and IVII ) . Actually, the EOS for this recipe 
has the unappealing feature that can be divergent in several cases [40j as we will show in Sec. IVII (c.f. Fig. |2T|) . This 
divergence is due to the fact that in a FRW cosmology the energy-density associated with this EMT becomes zero at 
some redshift (c.f. Figs. l22l and |23|) . a feature that had already been remarked by Starobinsky (25j |. 
Using Eq. ([5]) one can also write 



.4 



Kb ( A ) = -r /flflV a v fc i? 



fRRR{V a R)(V b R) - ?f ( Rf R + j 2/.T 



nT ab 



1 _ la 

A 



(10) 



Recipe III: The EMT defined in this recipe arises from Eq. ([3]) or Eq. ([5]) by simply identifying 



^ T T' X ■= fRR^aVbR + fRR R (VaR)(VbR) - 9ab 



~ Wr -f) + IrrUR + f RRR (VR) 2 



= f R RV a V b R + fRRR(V a R)(y b R) - 9 -f {Rf R + f + 2^r) . (11) 

so that Eqs. © and © read G ab = -fa ( T 1 ^ 1 ' x + T ab ^j . This EMT was considered by Sotiriou & Faraoni [H| 

(denoted T^fJ^ by them). It has the unpleasant feature that is not conserved, not even in the absence of matter, by 
the fact that Z^ 1 is not included as a factor in the r.h.s of Eq. (fTTI) . 

-,1V, X _ t-lrrllliX 



Recipe IV: In this case the EMT is defined by T ab ' := f R T l ab 11 ^ so that Eqs. © and © read G ab = 

b 



K ^ab X + K .f R l T a b- Like T^ b n ' , this EMT does not conserve either since in the presence of matter VT a ff x 



T ab ^ a \f R l ) ^ 0, however, unlike T I ' X , it is conserved in the absence of matter, as can be seen from the previous 



equation. As we shall discuss in Sec. IV Al several authors have obtained an EOS from the fourth recipe when applied 



to cosmology. In such scenario the previous equation will imply that the conservation equation associated with the 
geometric dark energy will contain a source term given by —T tt 'S7 t (f R 1 ) — pd^^^/dt (where t is the cosmic time) 
and thus such term will depend on the total matter density (i.e. radiation plus baryon plus dark matter densities). 
Incidentally, even if T ab ' ^ T b 1 ' X the EOS associated with both EMT will coincide as the factor f R x cancels 
when taking the ratio of pressure over energy-density for the X-component. 
When using Eq. © , T^ b v ' x can be written as follows 

1 A ' [fR R V a V b R + f RRR (V a R)(V b R) - (RfR + f + 2kT^ 



T J 

ab — „ f 

KJR 



(12) 



In our opinion Recipe I is the simplest, the best motivated EMT and the one that has the nicest features. In 
Sec. IV Al we shall derive the different EOS that arise from Recipes I-IV. 



IV. f(R) MODELS 

We have selected three f(R) models which have been analyzed carefully in recent years. First we shall present the 
models by reviewing some general features, and then test their cosmological viability using the approach that we will 
discuss in the next sections where we will compare our findings with previous results. In particular, we confront them 
in the light of the ACDM paradigm. 
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A. MJW model 



This model was proposed by Miranda et al. [301 ]: 

f(R)y U w = R-0R*^(^ + ^ S j , (13) 

where /3 and R* are free positive parameters, and it is well defined provided R/R* > — 1. In this work we use /3 = 2 
and R* = <t*Hq, where Hq is the Hubble constant today, and a* is a dimensionless parameter that can be adjusted in 
order to best fit the model to the cosmological observations. In the numerical analysis presented in Sec. I VII we take 
<7* = 1, which is not necessarily the best fit. Figures Q] and [5] depict the behavior of this model. The /—curvature 
Irr = PR* 1 /{R/R* + l) 2 (where we have dropped the label MJW) is positive definite, provided /3 > and R* > 0, 
and becomes small for large R (see Fig. [3J. In the domain where f(R) is defined, the derivative fa can be negative 
or zero in the range —1 < R/R* < j3 — 1, i.e., —1 < R/R* < 1 for the specific value (3 = 2. However, as we shall 
see in Sec. IVII it turns that fn>0 for the cosmological model that we analyzed, as the inequality R/Hfi > 1 largely 
holds during the cosmic evolution (c.f. Figs. [3] and |H] ), avoiding in this way any possible drawback associated with 
the values Jr < 0. Moreover, since R/Hq > 1, the pole at R/R* = —1, where fun — > oo and where f(R) is no longer 
defined, is never reached (c.f. Fig. 01, therefore, in this case it is unimportant to include Jrr in the potential V(R) 
that we introduce below. 

As we emphasized in the Introduction, the transformation from f{R) gravity to the corresponding STT can gen- 
erate several pathologies because the potential U(x) which is associated with the scalar field X — Ir mav be ill 
defined (notably when fun is not positive definite since then x(R) cannot be inverted in the original domain of R). 
Notwithstanding, the model proposed by Miranda et al. [3(| is free from those pathologies since in this case Jrr > 0, 
unlike the Starobinsky and Hu-Sawicky models discussed below. The potential V(R) introduced in Sec. [ill an d which 

is different from U{x) B, is given as follows for this model V{R) = ^ j(l + R){R + 6/3 - 1) - 2/3(3 + 2i?)ln(l +R)\, 

where R = R/R*. Figure [5] depicts this potential where the critical points (maxima or minima) at R\ correspond 
to a possible trivial solution R = R\ in vacuum (i.e. de Sitter point associated with i?i > 0), and are also the 
points, notably the minimum, that should be reached asymptotically, (in time for cosmology and in space for compact 
objects) when a non trivial solution for R approaches its asymptotic value in spacetime regions where the matter 
becomes practically absent or very diluted. The global minimum at Ri ss 6.1AHq is the actual de Sitter point reached 
in the cosmic evolution (c.f. Fig. [5]). Similar considerations regarding the critical points will also apply for the two 
models discussed below. 

Miranda et al. [3(| showed that such a model is consistent with a FRW cosmology and that singularities do not 
appear in idealized compact objects. Regarding these latter, we confirmed their findings using our approach [28j . As 
concerns the cosmological part, they specifically showed that their model is able to reproduce the matter-dominated 
epoch as well as the accelerated phase that are required to explain several observations. This is something that we 
confirm here (see Sec. lVI[) . This model was previously criticized on the grounds that, at the perturbation level, it seems 
to be unable to reproduce the observed matter power spectrum and that it may also be inconsistent with the local 
gravity constraints imposed by the Solar System experiments 42]. Miranda et al. 43] replied that the arguments 
put forward in Ref. [42j to discard the model are not strong enough and that a careful reexamination is needed. 



Furthermore, Thongkool et al. [44[ have also argued that this model is unable to satisfy the thin shell condition that is 
needed to produce a successful chameleon mechanism, being this presumably the only way to satisfy the Solar System 
constraints in f(R) theories. 

B. Starobinsky model 

The model is defined by the following function [25j |. 



f(R) s = R + XR, 



R 2 Y q 

1 + ^r -1 



R 



(14) 



4 The relationship between both potentials is given by dU(R(x))/dx = (2/ - f R R) /3 = dV(R)/dR, where U(x) ■= U(R(x))- Thus, 
dU(R)/dR = fjujdV(R)/dR, where we used d\/dR = fun- If we include the term fun in the definition of an alternative potential 
V(R) = (2/ - f R R) /(3/rr), then the relationship between U(R{\)) and V(R) dU(R{\))/dR = (f RR ) 2 dV(R)/dR. 
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FIG. 1. f(R) models considered here (see the main text for the analytical expressions). For reference the general relativity (GR) 
case /(-R)gr ~ R is also depicted as well as a GRA model /(-R)gra = R — 2A with A = 2.08Hq. Notice that for sufficiently 
large R/Hq, the Starobinsky and the Hu-Sawicky models can be approximated by their respective /(-R)approx := R — 2AJg. 
Such models behave as /(-R)gra (i.e. almost straight lines) with an effective positive cosmological constant given by the 

Value -/(0)approx/2. 




FIG. 2. Same as Fig. Q] for smaller R/Hg. 



with q, A positive parameters and Rs is another parameter playing the same roll as i?* in the previous model; we 
assume Rs — asH^, and take as ~ 4.17, q — 2 and A = 1 as in Ref. [37j. Both models, this and the previous one, have 
the property that, in vacuum, admit solutions with R = 0, like GR in vacuum (which includes the asymptotically flat 
solutions), unlike some models which contain a term ~ Rr 1 in the Lagrangian. It is by virtue of this property that 
Starobinsky entitled his model in terms of a "disappearing cosmological constant in f(R) gravity" . On the other hand, 
in the high curvature regime where \R\ 3> Rs, the model yields f(R) ~ R — XRs, an d thus it acquires an effective 
cosmological constant := XR s /2 (see Figs. Q] and |2J| , which is non-negligible provided A > 1. In this model, 
Irr — at R — ±R s /\/2q + 1, thus, fun is not positive definite (see Fig. 0}. The fact that this quantity appears 
in the denominator in Eq. Q indicates that a careful examination of the equations is required at those particular 
points, notably the positive one. Starobinsky himself ,25!], who dubbed such points weak singularities, has stressed 
the need of a close analysis of the solutions there. Nevertheless, we have not reached any of those weak singularities 
during the cosmological evolution as R > Rs/V§ ~ 1.86Hq (c.f. Figs. |4] and [9j. Figure |6] shows the potential 

1 /„, , _ 4R* + 5R 2 Ri + 3R%\ XRi 



] ' ( U ) = ^ ^ 2 - XRR S {R 2 +R S 2y ) + -^arctan(Ji/Jfc I '15, 
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FIG. 3. The derivative fu := df/dR for the f(R) models depicted in Fig. [TJ This quantity turns to be positive during the 
cosmic evolution. 




-1 ' 
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FIG. 4. The second derivative Jrr := d 2 f/dR 2 for the f(R) models depicted in Fig.Q] This quantity is positive definite in the 
MJW model and turns to be positive during the cosmic evolution for the Starobinsky and Hu-Sawicky models considered here. 



For R > 0, there is a global minimum at R = 0, a maximum at R ps 4.17iJg and a local minimum at R\ s» 6.82i7g 
which corresponds to the actual dc Sitter point where the cosmological solution settles in future cosmic time (c.f. 
Fig. [9]). Like in the previous model, the cosmological solution in the range we explored is such that fun and fu are 
always positive. It is somehow remarkable that AJ^ sa 2.09-ZJq is close to the actual effective cosmological constant 
A c g = R\/A sa 1.71-ffQ, even though the minimum is reached only for R\ sa 1.67i?s which one would think is not yet 
in the regime R^> R$- 

For the model with q = 1 and A = 1.56, we found previously static and spherically symmetric solutions that 
represent idealized compact objects with an asymptotically de Sitter behavior that was reached at the local minimum 
of V(R) (28j . In this regard, it is important to mention that we did not hit any of those weak singularities, notably 
the positive value) as the solution interpolates monotonically between R w 2.37i?sat the center of the object and 
R sa 1.98i?s asymptotically, which, as mentioned, corresponds to a de Sitter point However, it is fair to say that 
in that work we did not explore sufficiently the space of solutions and different values of the parameters. Moreover, 
we only used incompressible fluids. We plan to extend this analysis elsewhere. 

Starobinsky 's model can satisfy the conditions imposed by several cosmological observations |45|. For instance, 
it is able to produce an adequate matter epoch prior to the accelerated era, unlike several unsuccessful models (see 
Refs. [13, [U, HH for a thorough analysis). Moreover, the Solar System tests can be successfully passed by this 
model (e.g. taking n > 2) [25]. However it leads to ill defined potentials when transformed to a STT because fuR 
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X 

> 




R/H n 



FIG. 5. The potential V(R) associated with the MJW model Eq.([13]l with /3 = 2, R* = Hi. 



is not positive definite, and therefore the analyzes that rely on such an approach rise serious doubts about their 
soundness [24l. l26l. I27I I46l - l49l | . even if some of their conclusions turn to be true (c.f. Ref. [HI). 



Ed 
> 




R/Hn 



FIG. 6. Same as Fig. E] for the Starobinsky model Eq. (O with q = 2, A = 1 and R s w 4.17^ . 



C. Hu— Sawiky model 



The /(-R) for this model is given as follows pi 



f(R) =R-m 



2 Ci(E/m 2 )" 
c 2 {R/m 2 ) n + 1 ' 



(16) 



with n > and m 2 , ci and c 2 are parameters of the model. Like the previous two models, m 2 plays the same roll as 
i?* and Rs- According to HS1, the constant m 2 is fixed from the length scales of the Universe (it has the same units 
as R, so 1/m has units of length), 



77/ 



K 2 p 



(17) 
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where po is the average density of the Universe today. The numerical value that we assume and which is similar to 
the value taken in HS1 is m? ps 0.24Hq. The constants c\ and C2 are dimensionless parameters that can be fixed by 
demanding that this model mimics as close as possible the ACDM scenario. In particular their values are chosen for 
the model to match the current values ACCjU ^bar+DM « 0.24 and fi^v ~ 0.76 (see Sections fVl and |VT1 for definitions), 
where the upper-left index ACDM is used to distinguish from the corresponding matter content predicted by f(R) 
theories, which in general will not be exactly the same. The relative density will be replaced by a suitably defined 
energy-density of the geometric dark energy. Following HS1 one fixes such constants from 



ci 

(>2 



KCDMQO 

5 'baH 



DM 



and 



f°R - 1 = -n 



Cl 
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ACDMCtO 

"bar+DM 



- 9 



(18) 



(19) 



In this paper we take n — 4 and = 0.99 so that the specific values for c\ and C2 are c\ ps 1.25 x 10~ 3 , 
C2 fa 6.56 x 10 -5 . Figures [1] and [5] depict this model and show that in the regime R 3> m 2 , f(R) ~ R — cim 2 /c2, 
where now A 1 ^ — c\m 2 / (2C2) is not negligible provided c\ 3> c%. Like in the Starobinsky case, the Hu-Sawicky model 
has the property that fn and fun are positive during the cosmic evolution (c.f. Figs. [3l l4l and IT0| . even though they 
are not positive definite. 

The potential V(R) is shown in Figure[7] The analytical expression for V(R) in this case is not very enlightening (it 
is given in terms of a hypergeometric function) and thus we do not write it explicitly here, but as one can appreciate, 
its structure is similar to the potential of the Starobinsky model, except that the global minimum at R\ ps 8.9Hq 
corresponds to the actual de Sitter point found in the cosmic evolution (c.f. Fig. [TO]) . It is worth stressing that 
A^g ps 2.29H 2 while A e g = Ri/4 s» 2.23Hq. Those values are very close to each other since in this case R\ ps 37.08m 2 , 
which is already in the regime R 3> m 2 . 

As shown explicitly in HS1, the authors constructed a spherically symmetric static solution that represents the 
spacetime outside the Sun which passes the Solar System tests via a chameleon mechanism, as it provides a post- 
Newtonian parameter 7 pa 1. Furthermore, like the Starobinsky model, it is also consistent with the required matter 
dominated epoch prior to the accelerated expansion era as we shall see in Section IVII Notice that the Hu-Sawicky 
model with n = 2 is essentially the same as the Starobinsky model with q — 1 modulo a redefinition of their parameters. 

The Hu-Sawicky and Starobinsky f(R) models are perhaps the most tested models so far, which includes con- 
frontation with CBR, SNIa, baryon acoustic oscillations and gravitational lensing among others [UIHl 0- 



x 
> 




R/H n 



FIG. 7. Same as Fig.EDfor the Hu-Sawicky model Eq. (Q2). Here n = 4, m 2 = 0.24H-0 2 , ci w 1.25 x 10 -3 , and c 2 » 6.56x 10" 



5 The reader is urged to remember our definition of f(R) which contains the Ricci scalar unlike the HS1 convention. Therefore = 
1 + HS /fl- So our value /°, = 0.99 corresponds to HS /^ = -0.01. 

6 In several of theses references, generic f(R) models are tested using a parametrization that characterizes the deviations relative to GR. 
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V. COSMOLOGY IN f(R) 



We focus now on homogeneous and isotropic space-times which are relevant for cosmology and which are described 
by the FRW metric, 



where k — ±1, 0. 

From Eq. (|4]) we find 



ds 2 = -dt 2 + a 2 (t) 



dr 2 



1 — kr 2 



9 2 + sin 2 6 dip 2 ) 



(20) 



R = -3HR - 



3/ 



RR. 



SJrrrR 2 + V- IrR + 



(21) 



where ' = d/dt. Like in GR, the diagonal spacetime components of Eq. ([5]) will provide the remaining field equations. 
However, note that the t — t component of Eq. (JS|) will contain a term like R. This term will be rewritten in terms of 
lower order derivatives using Eq. (|2ip . The final expression for the equations which have the desired form of a Cauchy 
initial- value problem, subject to the initial data constraint is: 



H + — + T 

a 2 Jr. 



f rr.HR 



g VrR - f) 



3/fl 



(22) 



H 



where 



H 2 + ±- (f BR HR 

Jr. V 6 3 



H = a/a , 



(23) 



(24) 



is the Hubble expansion. 

Equation (|22|) is the modified Hamiltonian constraint which generalizes the usual Friedmann equation of GR. This 
equation constrains the possible values of a, H, R, R at some initial time. The explicit appearance of the scale factor a 
will come from the contribution of the energy-momentum tensor of the matter. We call these values the initial data, 
although this does not mean that it is the data at or near the big bang. Clearly, the above equations reduce to the 
standard equations of GR for f(R) = R. Moreover, taking f(R) — R — 2A one recovers the equations of GR endowed 
with a cosmological constant. 

Now, the expression for the Ricci computed directly from the metric is given by 



R = 6 [H + 2H Z 



(25) 



Note that by using Eqs. (|2"21 and (|2"B")) in Eq. (|2"5)) we obtain an identity R = R, which shows the consistency of the 
equations. Thus, as emphasized before, this equation is redundant. However, one has the freedom of using the system 
of Eqs. (|2Tj). ([22]l and ([25]) instead of Eqs. (|2T ]l -([2"3" l) to solve for (H,R), and Eq. {24]) to solve for a. In our case, we 
have in fact, used both sets of equations in order to check the consistency of our computer codes and also the accuracy 
of the numerical results at every "time" step (see below). We stress that Eq. (|2"2"|) is only used to fix the initial data 
and then it is monitored to check that it is fulfilled at every integration step within the accuracy of the numerical 
algorithm (4th order Runge-Kutta scheme) when integrating the dynamical equations. At this regard, Motohashi et 
al. (3(| had also used a similar technique to check the numerical accuracy. 

As we mentioned before, the matter variables T a ^ obey their own dynamics which is provided by V ' a T ab = 0. We 
shall assume that T a t, is a mixture of three kinds of perfect fluids, T a t, — Yli=i ^ab> ^ab = (P» +Pi) u aUb + gabPi, which 
correspond to baryons, radiation, and dark matter, respectively, in a epoch where they do not interact with each 
other except gravitationally. The energy-momentum tensor for each matter component is conserved separately, and 
the conservation equation V a T i ab = leads to the usual expression 



-3H(pi +p 4 



(26) 
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where the total energy-density of matter is p = Y^i=i Pi- The above equation integrates straightforwardly (with 

Pbar,DM = and p rad = Prad/3) as follows 



Pbar + PT)M , Prad 



(a/a ) 3 (a/ao) 4 



(27) 



where the knotted quantities indicate their values today. Thus in this case the matter variables that appear in the 
field equations are given explicitly by T\ = -p and T = £) i=1 T. t = £) i=1 (3p» _ = _ (Pbar + Pdm) (since T rad = 0, 
and Pbar.DM = 0), which in turn can be written in terms of the scale factor according to Eq. (I27|) . In this way, the 
differential equations will depend explicitly on a(t). 

In the majority of the cosmological studies performed so far the equation for H is usually written in terms of R 
which is clearly not suitable for an initial value problem (e.g. see Ref. [53| for a system of equations of this sort used 
for reconstructing f(R) functions), unless one uses Eq. (|2"TT) to replace such a term in favor or lower derivatives or 
alternatively if one uses a combination of several variables (e.g. see Sec. 4.1 of Ref. [l3j], and [H|). Our approach is 
similar to that of Appleby & Batty [54| where they solve the set of Eqs. (|2"Tj) . ([24]) and (|25|) . but it is different from 
theirs in that they do not use Eq. (|22[) to fix the initial data and monitor the consistency of the numerical solutions 
(see Sec. IVBp . In Sec. IV Bl we shall rewrite the system of equations using another "time coordinate" which is more 
suitable for the numerical integration. This is another important difference with respect to Ref. [lij where the authors 
use t itself as independent variable and furthermore they integrate backwards in time, which as we argue in Sec. IV Bl 
may have several inconveniences. 



A. The Equation of State in f(R) 



Let us now study the specific application of the EMT's considered in Sec. IIIII In GR with matter and dark energy 
components, the parameter a is directly related to the total EOS cj to t = Ptot/ptot [i.e. a/a = — («;ptot/6)(l + 3w to t)], 
which in turn, determines if the universe is expanding in accelerating or decelerating way, provided ptot > 0. Here 
Ptot = P + Pde, and ptot = Prad + Pbe- One can then define wde = Pde/pde, as the EOS for the dark-energy, which 
will coincide with a; to t when the matter contribution is negligible as compared to that of the dark-energy; this occurs 
in the ACDM paradigm at future time. 

We can in a similar way define an EOS for the component associated with f(R) gravity. In order to achieve our 
goal we need then to define the energy-density px and pressure px, associated with the geometric dark energy fluid. 
We shall obtain three expressions corresponding to Recipes I-IV. 

Recipe I: First, we define the energy-density px so that the modified Friedmann Eq. (|2"21 reads (hereafter we 
asume k = 0) 

H 2 = ^(p + p x ) , (28) 

where we remind the reader that p includes all the contributions of ordinary matter (baryons and radiation) and dark 
matter as well. Second, we define its pressure px so that Eq. (|2"5)l reads 

H + H 2 = -^{p + Px +3(p rad +px)} • (29) 

In this way, Eqs. (|22|) and (|23|) together with the above two equations lead to the following expressions: 

Px = \\UrR - f) - RB.HR + «p(l - fn) \ , (30) 



Px = -^r- ( \ UrR + f)+ Sf RR HR -k( P - 3p rad / fl ) } . ( •'! I ) 



The EOS ux of the X-fluid reads then 



ojx = — ■ (32) 
Px 
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The above equation can also be written in the following implicit form when using Eqs. (|25|) . ([28]) and (|29|) 



OJX 



3H 2 - 3n Prad - R 
3 {3H 2 - up) 



(33) 



Clearly we have assumed here that f(R) ^ as otherwise one is led to ujx = 0/0. Numerical examples of this 
EOS will be provided in Sec. EJ for the three f(R) models introduced in Sec. [TV] (c.f. Figs. [T8H20J). 

One can recover Eqs. (JSOJ) and (JSU) directly from Eq. © using p x := u a u b T*, with u a = (d/dt) a ), p x ■= S x a a /3, 
where S x b is the 3-energy-momentum tensor obtained from and which is denned on the orthogonal hypersurfaces 
to u a El, and by replacing R (which will appear in px via the term VtVt-R) in favor of lower order derivatives using 
Eq. (EH). 

As we have emphasized, T x h is conserved V a T^ = 0. Indeed, since the only energy-momentum tensor compatible 
with the hypothesis of homogeneity and isotropy is that of an effective perfect fluid, then a fortiori px and px must 
obey an equation similar to Eq. (f2l))) . This statement can be checked explicitly by using Eqs. (|3"0|) and pip in Eq. (|26l) . 
and then appealing to the field equations. 

The other general considerations discussed in Recipe I shall apply to this particular case. For instance, although the 
matter terms appear in Eqs. (1501) and (j3"Tj) . px = = Px for the GR case f(R) = R, and moreover, px = A/re = — px , 
lux = — 1, for f(R) = R — 2 A. This shows that these quantities so defined are sensible because it is only when the 
theory differs from GR that they do not vanish in general, and when A is included, they reduce to the expected values. 
Of course, we do not intend to add a cosmological constant in the f(R) models, the only purpose of these remarks is 
to illustrate some of the properties of these definitions. Notice that these expressions do not correspond to the total 
energy-density and pressure of the effective "matter", i.e. the sum of ordinary matter, dark matter and A-matter, 
but only to the A-fluid. The total equation of state is given by 



Ptot 
Ptot 



(34) 



where p to t = Px + P [c.f. Eq. (gll)] El and pto 



Px +Prad [c.f. Eq. (|29|)], The total EOS is depicted in Fig.lH These 
quantities can also be computed directly from T*£* as ptot — u a u b T^ and p to t ■— T t ° t ' t /3. 

As part of Recipe I, one can write px and px in terms of "purely geometric" quantities where the matter terms do 
not appear explicitly. This is in fact achieved by considering the equivalent definition Eq. ([7]). Let us see explicitly 
how can we arrive, for instance, to px directly from the field equation, which is easier than using Eq. ([7]). Consider 
Eq. (|22[) , with k — for simplicity, and again assume a perfect fluid. Then we have 



1 



= -3H 2 f R + - (f R R - /) - 3f RR HR + k P . 



(35) 



Adding 3H 2 to both sides we then write 



H 2 = 



3H 2 (1 -f R ) + l (f R R -!)- 3f RR HR 



(36) 



Therefore, this equation reads like Eq. (|28|) if one defines the term in brackets as follows 

npx - 3ff 2 (l -f R ) + l (f R R - /) - 3f RR HR . 



(37) 



In this way, the matter term does not appear explicitly in px, but clearly (1371) and (|30[) are one and the same thing. 
We can recover the previous equation if one uses up = 3H 2 — npx in Eq. Q30[) and then solves for px- One can rewrite 
px in a similar fashion (c.f. Ref. (30]|). 

Miranda et al. [3(| , de Felice & Tsujikawa [l3[ , Motohashi et al. [36l - l38l | , and Bamba et al. [55| have used a definition 
for the EOS which coincides with Recipe I. In the case of HS1 [llj], their definition coincides exactly with our expression 
Eq. (|33p when the radiation contribution is completely neglectecF^I. which is justified only in the matter dominated 



7 Of course f(R) = C = const, is also excluded as otherwise one is lead to consider a vanilla gravity described by Cg a t,/2 = —kT^. 
Moreover, this would lead to fn = and px and px lead also to 0/0. 

8 S x a a = i is the spatial trace of T\. 

9 Under this definition f! := Kp to t/3H% = 1 for k = [c.f. Eq. (57)]. 
In HS1 the authors use completely different variables. 
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epoch. It is to be notice, that some authors have preferred to keep the terms like R and a in their definition, instead 
of writing them using lower order derivatives as we do. 

Equations (j3"0)) ~(l3"2" j) or even Eq. (|3"3"|) . are rather simple and easier to calculate under our approach. They are 
evaluated easily at every integration step when solving numerically our system of equations (see Sec. IVB"|) . 

Recipe II: In order to obtain p^ associated with Recipe II we can use Eq. © or (fTU|) . However, for our purposes 
Eq. (fTU)) is more convenient. We obtain then 



Px=^r{~ UrR - f) - ZIrrHR + kp I 1 - ^ 



A 



A 



(38) 



n 11 

Px 



2 (IrR 



/) + ZIrrHR 



. > f R 
oPrad^j- 



(39) 



Equation (I3"51) can be recovered easily from Eq. 



we define p 1 ^ from F*l 



following almost the same steps as in Recipe I, except that now 



AH 2 



(p + pT 



(41) 



Therefore the EOS reads 



n ii 
ii _ Px 
u x - — 
Px 



(42) 



The quantities of Recipe I are recovered for A = 1. In Refs. [2l|, I39l441j Recipe II was used with A = f R (F in their 
notation; see Ref. [4l[ for a discussion concerning the introduction of this constant /£). 

In general f R =^ 1 although f R ss 1 (c.f. Ref. [ill) H As shown by Hu & Sawicky fig , one has to take (at least in 
their model) ss 1 to avoid large deviations in the EOS relative to the cosmological constant value uj\ = — 1. So if 
one takes f R ^ 1, wj^ and wx are not equivalent. As a matter of fact, wj^ can diverge at some redshift depending 
on the f(R) model. For instance, in the Starobinsky and MJW model this is exactly what happens (c.f. Figure [2~Tj) . 
because p 1 ^ becomes zero (c.f. Figs. 1221 and |23|) . In the Hu-Sawicky model we did not encounter that divergence in 
the range of redshifts we explored (c.f. Fig. [20]) . Due to this pathological behavior, this definition for the EOS makes 
it rather unsuitable for comparing with observations. 



Recipe III: In this case p^ 1 and p^ 1 are obtained from Eq. (TTTj) . In particular, p 1 ^ 1 can be read off directly from 
Eq. m, by defining H 2 = + p)/(3f R ), and get 



PW 



(f R R-f)-3f RR HR 



(43) 



The pressure and the EOS read respectively 
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(44) 



and 



in _ Px 
,J x — —m i 
Px 



(45) 



11 Eq. I|38| l has the alternative expression 

K P ][ = 3H 2 (A - / K ) + ~ (f R R - /) - 3f RR HR . (40) 

The interested reader can consult Ref. [4l| for the "purely geometric" expression for p 1 ^ . 

12 Since this dimensionless quantity appears as G e g = Gq/ fu and at present time Gg ff ~ Go, thus one expects f ( R Ri 1. 
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where T = 3p ra d — p, as before. This EOS coincides with lux only in vacuum. But even in vacuum the density and 
pressure does not satisfy a conservation equation like (p?o| . as it was stressed in Sec. IIII1 



Recipe IV: In this case p 1 ^ and p 1 ^ arise from Eq. (TT2)) or from the definition T , ' :— /p 1 T ' I ' as given in 
Sec. IIII1 which yields 

(46) 

JR 
„III 

PW = ^- (47) 

JR 

The EOS is thus 

uj x - —jy = lu x . (48) 
Px 

Alternatively, p 1 ^ can be read off directly from Eq. (|22jl. by defining H 2 = k (p™ + p/ Jr) /3. 

Notice that the Recipe IV quantities coincide with px, px and lux only in vacuum. Nonetheless, in the non vacuum 
case p l x and p l x have the unaesthetic feature of not satisfying a conservation equation like f|26[) , but rather present 
an extra source term that depends on the total matter density r 3 ]. Quantities associated with Recipe IV have been 
considered by Sotiriou & Faraoni [l2j E3, Capozziello et al. [HI l56l| l 15 l. and Bamba et al. [E3] , although their expression 
for the pressure is written in a different way (a term R appears explicitly there). 

Despite the unpleasant features associated with lu^J 1 — (j™ \ they can behave similarly to lux for actual cosmological 
scenarios (c.f. Figs. fT8HT9|) . 

In summary, there exist at least three inequivalent definitions lux, wj/, 10 x 1 that we have identified in the literature 
as arising from different sorts of EMT representing the geometric dark energy component. We consider that lux has 
the most appealing properties, while EOS ou x may have serious deficiencies as it can diverge at some redshifts. On 
the other hand, wj/ 7 = lu x have the unaesthetic feature of being related to a non conserved EMT. But even if we 
dismiss this feature, there are quantitative differences between lux and lu 1 x 1 (or ui x ) ano - although the numerical 
discrepancies between them maybe considered unimportant at present time for the viable cosmological models, they 
may, however, become important if the EOS is measured with high precision in the future. In this latter instance, 
one should bare in mind which definition is being used to make comparisons with the values that are inferred from 
observations. We cannot emphasize more this last comment. 

In the following sections we describe the numerical strategy we use to integrate the system of equations (indepen- 
dently of the specific f(R) model considered) and then analyze the cosmological results for the three specific models 
described in Sec. IIVI 



B. Numerical integration 

The system of Eqs. (|2l]). ([23]) and (J24]) as well as the alternative system (|24|) and (JUJ have the form dy i /dt = 
F l (y % ) where y l — (a, H, i?,II) and II := 7?. Therefore they can be solved easily with a fourth order Runge-Kutta 
algorithm. Now, it is well known in cosmology that t is not the best independent variable to perform the cosmic 
evolution because usually several scalars blow up very fast as a — > 0. It turns better to use the following parameter 
to integrate the differential equations F^l 

a = ln(a/a ) , (49) 



3 The explicit form of the source term can be appreciated in Eq. (108) of Ref. [TJ], Eq. (8) of Ref. [5(| and Eq. (11.3) of Ref. [5^1 . 

4 In Sotiriou Sz Faraoni [T3 | it was defined an EMT denoted as T^J ^ which corresponds to our T*l I,X . Incidentally, those authors 
define also p e ff and p e ff which correspond to our p 1 ^ and p 1 ^ , rather to p 1 ^ 1 and Pjf 1 ■ This can be appreciated in Ref. ll'J as the 
quantities p e 1 1 and p e 1 1 do contain the the alluded factor , which is missing in their T^l ff) . It seems that those authors defined 

p e ff and p e f f using the cosmological field equations instead of using T^J^^' directly. In any case, this apparently lack of consistency in 
their notation has no effect on the EOS as Recipes III and IV give rise to the same EOS. However, the reader is urged to have in mind 
such nuances in order to avoid any confusion between the current paper and theirs. Still, none of the X-fluid quantities associated with 
Tab 1 ' X or '^ab' ' w '" satisfy a conservation equation like Eq. J26D . 

5 A typographical error (which seems to be systematic) is found in Capozziello et al. Refs. [Til. |56H where a sign ' — ' seems to be missing 

in front of the term i (/ — }rR) when they define p CU rv = -^j-^ (/ — fnR) — SfuRHR^ [c.f. Eq. 1 143 I I and Eq. (1221 in vacuum 

p = = -T\ ] . This is not a global sign, so their definition for p C urv seems to be inconsistent with their modified Friedmann equation 
which (modulo notation) should be identical to ours. Furthermore, taking f(R) = R — 2A yields p C urv = —A/k which has the opposite 
sign that one should usually consider for A > 0. The expression foj j) curv is of course invariant under / — > — /, thus the difference cannot 
be explained by a sign convention in the gravitational Lagrangian. 

6 Clearly, a is a good "time coordinate" provided a(t) is a monotonic function, otherwise the relationship between t and a is not in 
one-to-one correspondence, making the mapping t(a) ill defined. 



which maps the big bang a — > to a — > — oo. So when integrating the equations with respect to this variable one 
is always far from the big bang in the a domain, but one can be very close to it in the t domain, in the sense that 
a/ ao <C 1. Equations (|2"Tj) , (f2"2")l , (j2"3")l . and ([25]) read as follows with respect to a: 



R" = -R' ( 1 



R 



[3fRRRH 2 R' 2 + 2/ - f R R + kT] , 



(50) 



(51) 
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fR R H 2 R , --(f R R-f) 



3/fl 



(52) 
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JrrH 2 R' 
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G 



(53) 



As we have mentioned, we can use Eq. (IBlj) or Eq. (|53|) to solve for H. Nevertheless, we have used both equations to 
verify the consistency of our numerical code and together with the verification of the modified Hamiltonian constraint 
at every integration step, prove the soundness of our results. 

From (gj), and one obtains 



dt 



(54) 



which is used to extract the age of the Universe when inverting the numerical solution t(a) so as to obtain a(t) 
(the age is "read off" from its graph). The initial condition needed to solve this equation simply sets some irrelevant 
initial time which can have any arbitrary value. The relevant chronological quantities are the differences between any 
given time and the initial time. 

Different system of equations us ing different variables have been analyzed in the past in order to perform the 
integration (e.g. see [H, [H, [2l|, [58| ) . In particular, the system employed by HS1 [l8[ has become very popular in 
recent years. 

We stress that the numerical integration must be performed forward in a because there are solutions in f{R) that 
behave like an attractors as one evolves towards the future. Integrating in the opposite direction can then easily 
lead to cosmological solutions that do not satisfy (retrodict) the physical conditions they should have in the past 
(e.g. those dictated by CBR or nucleosynthesis). We think that this is precisely the problem that was encountered in 
Ref. [54^ . where some kind of "singularity" was found due to runaway solutions which corresponded to "inadequate" 
initial data. Such data, that we can call dfj and which is fixed in the "future", might be slightly different from the 
data that is predicted when integrating forward in time using d° as initial data fixed in the past, but this difference 
is large enough not only to prevent the system from retrodicting the initial data dp, but also to make it incapable of 
arriving to the same past epoch associated with dp. One would then require a remarkable precision on to retrodict 
dp] such are the features of dynamical systems that have attractors in one direction of the evolution parameter. It is 
beyond the scope of the present paper to analyze such dynamical properties of the system. 

The numerical strategy we used to integrate the equations is as follows : We start the integration at some initial z 



in the past (z c ), where z = ag/a — 1 



1, given initial conditions for p" 



„rad 



R c , R' c , H c . From Eqs. (JSU) 



and (|33[) evaluated at the initial z c , one solves algebraically for R' c and R c , respectively, in terms of the remaining 
variables. In this way, both R c and R' c are fixed given the initial values H c , p" latt , p r c ad and i>j c x . We can fix p™ att 
and p T c ad using the ACDM as a guide to chose the abundances of baryons, DM and radiation at that z c , and take 
uj c x w — 1. This value for ui x might seem rather ad-hoc at first sight, however, since one expects that for "high" z the 
curvature R is also sufficiently high for the f(R) models to behave as GR plus A c g (c.f. FigfTJ then one can anticipate 
that lj c x w — 1 is not a bad approximation. Moreover, we have adopted such value in order to compare our results 
with other authors who have taken similar values for the EOS at the same z c . In any case, the value for lu x is not a 
fundamental issue and initial conditions can be fixed as one pleases provided the modified Hamiltonian constraint is 
satisfied and that the resulting cosmological model is consistent with observations within the error bars. So, what we 
have used as main criteria to asses if of our initial data is adequate is that our results at present time be consistent 
with what is measured from observations. Finally, the remaining quantity to be fixed for the rest of the initial data 
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to be determined is H c . The AC DM model of GR provides also a clue of the expansion rate at a given z, and in 
order to fix it more accurately we used a shooting-like method for H c that allowed to recover the actual abundances 
for each model at present time when integrating the system (|50p . (jSTj) and (|55| . As we have stressed, Eq. serves 
only as a constraint for the initial data and also to monitor the accuracy of the numerical results at every time step. 
This equation is extremely sensitive to any typing mistake in the computer code as well as to any inconsistency when 
fixing the initial data. A final comment is in order concerning the different recipes for the EOS. Since we use recipe I 
given by Eq. to help us fixing the initial data, we do not expect that the initial (and in general any) values for 
the other EOS, wj/, and wj/ 7 will be the same as u)x- Alternatively, we could have opted to fix the initial data so 
that the three EOS initially behave the most similarly, knowing that proceeding in this way could ruin the prediction 
for the present values of the observables. This method can certainly be followed, but we do not pursue it here. The 
important point to keep in mind is that our resulting cosmological models are viable given a self-consistent initial 
data. 

It is important to point out some differences about the fixin g of initial data with respect to some recent works that 
have appeared in the literature. For instance, several authors |56l. l59l Io0| have proposed to use the current values of 
the deceleration parameter and the jerk, respectively related to do and a'o, to fix the initial conditions for Rq and Rq. 
For instance, such quantities can be written respectively asFl 
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aH 2 
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iH 3 6H 3 



H 2 + H 
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H 2 H 
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' QH 2 

R 
6H 3 ' 



(55) 
(56) 



This would be a nice strategy if one could easily integrate from present to past, nevertheless, as we have emphasized, 
one easily goes into inadequate solutions by proceeding that way. Following our integration strategy, the deceleration 
and jerk parameters are rather predicted from our initial conditions imposed in the past. In the next section we 
provide today's values of these parameters for each of the three f(R) models analyzed here. 

VI. NUMERICAL RESULTS AND COSMIC VIABILITY 

In this section we present the numerical results for the three f(R) models introduced in Sec. III. Figures [8lI10l depict 
the Ricci scalar R as a function of the light-shift z = ao/a — 1, where z = corresponds to the present time. For the 
three models, R approaches the de Sitter minimum of the potentials plotted in Figures EHT] The expansion H damps 
R, and near the minimum, it oscillates. The oscillations can be easily explained by a linear perturbation around the 
de Sitter minimum where the matter contribution is negligible as it dilutes rapidly with the expansion. 
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FIG. 8. Ricci scalar for the MJW model as a function of the light shift z. 



Notice that some authors define the jerk with the opposite sign [6111 
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FIG. 9. Same as Fig. [8] for the Starobinsky model. 
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FIG. 10. Same as Fig. [8]for the Hu-Sawicky model. 

Figures [TTUT31 show the expansion rate H. Like in Refs. [Hj], we see that the Hubble expansion oscillates at late 
times for the three models. 

We introduce now the dimensionless densities of the different "species" as follows: 

fl := O rad + fibar + ^DM + = 1 , (57) 

where fli = Kpi/(3H 2 ) 0- Here px is taken as in Recipe I. 

Figures [TH - [TB1 depict the relative abundances Qj as a function of light-shift z. For reference, the corresponding 
abundances of the ACDM model are also plotted. The abundances in the ACDM model have the following an- 
alytical expression in terms of their values today (knotted quantities) and the scale factor a = a/a : Q ACDM = 

ACDM Q°a I [( ACDM ^ ar + ACM ^dm) a~ 3 + ACDM ^d^ 4 + ^aV\ wher e the subindex i stands for radiation, 
baryons, dark matter or dark energy (i.e. the cosmological constant), and the exponent, /, takes the following values 
/ = —4,-3,-3,0 for the previous species respectively. Notice that Xa=i ^ cdm = 1- The current abundances fl® 
associated with the f(R) models are predicted given the initial conditions in the past, since, as we emphasized above, 
we integrate from the past to the present time. Although the initial conditions may vary from model to model, we 



Some authors 39-41, 45] define f!; = Kpi/(3fRH 2 ). This difference must be bare in mind when comparing results. 
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FIG. 11. Hubble expansion for the MJW model as a function of the light shift z. 
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FIG. 12. Same as Fig. QT]for the Starobinsky model. 

try to fix them in order to predict the actual abundances today. Nevertheless the fixing is not exactly the same for all 
the models, and so the current abundances will change slightly from model to model. We have taken into account the 
radiation contribution, however, it is almost negligible during the baryon-DM and dark energy domination eras. We 
notice that the Starobinsky and Hu-Sawicky models behave very much like the ACDM model of GR. In particular, 
the Hu-Sawicky model is almost indistinguishable from the ACDM model. However, the MJW model shows some 
differences at higher redshifts. 

The three models exhibit an adequate matter domination era, at least in the range of redshifts explored in the 
numerical evolution, which is followed by an appropriate accelerated expansion afterwards. This is a very important 
test in view of the results of Amendola et al. [20|, |2lj, l23[ who showed that several f(R) models are simply unable 
to recover an adequate matter dominated epoch or a suitable accelerated expansion. Furthermore, we see that fix 
decreases like Q\ for large z, and if we extrapolate this behavior to the very early Universe, it is expected that 
radiation will dominate and hence, that the standard predictions of primordial nucleosynthesis will not be spoiled B 

Figure IT71 depicts the behavior of the scale factor during the cosmic evolution for the three f(R) models and it is 
compared with the ACDM model of GR. The curves show the initial epoch at which the equations started to be 





If one does not want to extrapolate any behavior whatsoever, then one needs to integrate the equations from the nucleosynthesis time, 
where radiation dominates, to the current time, by fixing the initial conditions at nucleosynthesis by demanding a successful abundance 
of the primordial light elements and in addition a successful matter and dark-energy domination at late times. This cannot be an easy 
task, but it may certainly be performed. 
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FIG. 14. Densities for the MJW model as a function of the light shift z. For reference the densities of the K.CDM model 
are also plotted. Radiation has been taken into account but cannot even be noticed during this epoch. Here we assumed 
ACCM ^tar+DM ~ 0.24, ACDM fi° ad w 4.1 x 1(T 5 and Q,° A sa 0.76 for the ACDM model of GR. The predicted values of the 
corresponding densities today for the f(R) model are slightly different as they are predicted given the initial conditions in the 
past. 



integrated (see the discussion at the end of Sec. IVI A[) . The figure shows that in the Starobinsky and Hu-Sawicky 
models the age of the Universe is of the order Hq 1 ps 9.78/i _1 x 10 9 y (where h = (i?o/100)km _1 sMpc). While in the 
case of the MJW model, the Universe is slightly younger, 0.94iJ _1 pa 9.19h~ 1 x 10 9 y. Taking h = 0.7 as in Figure l27l 
we obtain an age ~ 13.97 x 10 9 y for the Starobinsky and Hu-Sawicky models and ~ 13.13 x 10 9 y for the MJW model; 
both agree with estimates of globular clusters in the Milky Way [62j . 



Figures [T8r[2"0l show the behavior of the EOS associated with the X -fluid based upon the Recipes I— III. As concerns 
the Starobinsky and the Hu-Sawicky models, our results are consistent with those reported in Refs. [H, HB-S, [55| 
for similar values of the parameters where Recipe I was used. We also obtain similar results for the MJW model as 
compared with Ref. [3(|, where Recipe I was used as well. The inequivalent definitions of the EOS given by Recipes 
I and III give similar results for the three models. However, the EOS of Recipe II can have a completely different 
behavior as it diverges at some redshift for the Starobinsky and the MJW models (see Fig.l2ip. This divergence, which 



was reported in Ref. [40j and commented in Ref. [25[ , is due to the fact that px becomes null and then negative, as it 



is explicitly shown in Figures |2"21 and |2"51 This in turn can be understood by looking at Eq. (J38J, notably at the term 
p(l — In/ A). First, in the case where A = 1, which corresponds to px, the term p(l — fn) is always positive since 
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FIG. 15. Same as Fig. [T4]for the Starobinsky model. 
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FIG. 16. Same as Fig. [14] for the Hu-Sawicky model. 



Q < f R < 1 during the cosmic evolution (c.f. Fig. [3]). If the other contributions are positive, which is expected since 
Irr *C 1 and since the other terms will give rise to A e g > 0, then px > 0. Now if A = f R < 1 which corresponds to 
Px , the term p(l — fR/f R ) can be negative in epochs where f R > f R which are precisely those epochs of large R and 
large z. In the matter domination era the negative term p{\ — /r//r) can dominate over the geometric /(i?)-terms of 
Eq. (|38|) . This is exactly what happens as depicted in Figures [22l and [23l As the evolution continues, p decreases to a 
point where the other terms of Eq. (|38[) . which increase, balance exactly to give Px = 0. As the evolution continues 
to a point where f R < f R or p is sufficiently small, the term p(l — fR/f R ) becomes positive or small, and then p 1 ^ 
becomes positive. All this behavior exacerbates as the parameter f R differs sufficiently from unity, since then the 
term p(l — fR/f R ) becomes important. Since in the Starobinsky and MJW models, unlike the Hu-Sawicky one, the 
parameter f R is not fixed in advance but rather predicted from the initial conditions in the past, then for those models 
it turns that f R differs from unity more importantly than in the Hu-Sawicky case where the model was constructed 
in such a way that at present time f R « 1 (see Table HJ. In the Hu-Sawicky model then the term p{\ — fR/f R ) is not 
very important in the epochs where p is not very large since the factor 1 — fR/f R ~ 0. However, farther in the past 
where p dominates, the term p(l — fR/f R ) might become important and can make p^ = in that model too. This 
depends on the rates at which the factor 1 — fR/f R and the energy-density p decreases and increases respectively as 
looking from the present to the past. In any instance, this behavior is unacceptable and indicates that the definition 
of EOS II is not adequate. 

Like in Refs. [TH, I3614381 loll [55| . we also find oscillations of the EOS around the "phantom divide" ujx = — 1, for 
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FIG. 17. Scale factor as a function of cosmic time for the MJW, Starobinsky and Hu-Sawicky models. The ACDM model is 
included for reference. Here to = Hq 1 . The age of the Universe can be estimated from the time when a/ao = 1 (today) and 
the time when a/ao ~ ("big bang"). 
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TABLE I. Deceleration, jerk and other quantities of the models at z = (today). 



the Starobinsky and Hu-Sawicky models. However, for the MJW model, the well behaved EOS, ux and do not 
cross the phantom divide in the past. This behavior was already remarked in [30j . 

Figure depicts the total EOS (i.e. using Recipe I). We appreciate the transition from the matter dominated era 
to the accelerated era. Figures l2"5l and 121)1 depict the deceleration parameter and the jerk as computed from Eqs. ([55)1 
and (|56p for the three f(R) models. Table U shows the current values of these and other quantities for such models. 
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FIG. 18. Equations of state according to Recipes I — III defined in the main text for the MJW model. 



24 



-0 . 8 



-1 . 1 



-1 . 2 





\ 


co x (I) 

co x (II) 

co x (III) 


1 














/ 




\ 




/ 

/ 






k / 

\ v 





FIG. 19. Same as Fig. [T8]for the Starobinsky model. 



-0 .99 



-1 .01 



-1 .02 







(O x (I) 


/ 

/ 


X 


C0 X (II 


j- ; 

I) 


J 


\\ 

\\ 








Va 

\ 



























FIG. 20. Same as Fig. [18] for the Hu-Sawicky model. 

A. Luminosity distance and SNIa 

The accelerated expansion of the Universe was corroborated by the measurements of luminous distance in supernovae 
(SNIa), when comparing with the predictions of the AC DM of GR. In which follows we confront the same data with 
the accelerated expansion predicted by f(R) theory. We compute first the luminous distance given by 



.at = m 



where a = a/ao and 



( = cH~ 



da* 



a a* 2 iJ(a*) 



(58) 



(59) 



where we have introduced explicitly the speed of light c in order to compute the distances in units of Mpc, and 
H := H/Hq. Another useful quantity in cosmology is the angular diameter distance given by 

DT = aC(a) . (60) 

We emphasize that the above expressions for the luminous and angular diameter distances are valid only for k = 0. 
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FIG. 21. Equations of state according to Recipe II defined in the main text for the Starobinsky and MJW models. Notice 
the divergence due to the corresponding energy-density becoming zero (see Figs. 1221 and [ 23)) . For lower z these EOS are also 
plotted in Figs. [T51 and [191 
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FIG. 22. Energy-density px (solid line) and pressure Px (dashed line) computed from the MJW model using Recipe II. The 
quantities are given in units of the today's critical energy-density. Notice that px becomes zero at z ~ 2.38. 



Perhaps the easiest way to compute £ within the framework of the numerical scheme that we have devised is to 
transform Eq. (|5^|) into the following differential equation for 



da 



1 



a 2 H{a) ' 



(61) 



where £ = £/(c H 1 ) is dimensionless and which in terms of the variable a = ln(a) defined in Eq. 



reads 



C' 



H 



(62) 



This first order differential equation is integrated simultaneously with the field equations in the way we just described. 
A technical but important point is that a priori we do not know what the initial value of £ is at the past epoch where 
the numerical integration starts. However, that value is easily found by a shooting method such that £ is zero today 
(at z = = a). Given £ one then computes d^ at and Z?^ at from Eqs. (|55|) and (|6Tj|) . Actually the quantity which is 
usually reported is not d^ at but the distance modulus given by 



fj, 



M = 51og 10 (d« at /Mpc) 



(63) 
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FIG. 23. Same as Fig. [22] for the Starobinsky model. In this case p x becomes zero at z ~ 2.1 
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FIG. 24. Equation of state cjtot as defined in Eq. 



Figures l27l and l28l show the luminosity distance for the three f(R) models described above and are compared with 
the ACDM model. Here the data corresponds to both the "historical" data of Riess et al. of Ref. 0] and also the 
UNION 2 data of Amanullah et al. Q . We appreciate that up to z — 1 there is no significant difference between the 
three f(R) models and the AC DM model. Nonetheless, for larger z (see Fig. [29)) . significant differences arise which in 
the future could constrain or rule out these or other f(R) models or even challenge the ACDM paradigm. Figure l30l 
depicts the angular diameter distance. 

To conclude this section we mention a technical point that in principle restricts the integration of the equations for 
some f(R) models and which manifest specifically in the Starobinsky and Hu-Sawicky case where fuR approaches 
zero for large R/Hq. This regime occurs during the cosmic evolution for large z. Figures [31] and [32] depict the behavior 
of fim, where we appreciate that frtrtHQ <C 1, and approaches zero much faster in the Starobinsky and Hu-Sawicky 
models than in the MJW model. This behavior is consistent with Figures Q] and 0] where one sees that for large R the 
Starobinsky and Hu-Sawicky models behave as f(R) ~ R—2A e g, and therefore frtR ~ 0. Now, the point is that when 
fnR ~ in the past, Eq. ([50l) develops large variations as f RR appears in the r.h.s and affects the numerical precision. 
Since in the MJW model fRR decreases slower than in the other two models, we have corroborated that for the former 
we can perform the integration starting quite far in the past (even at the recombination epoch) without finding any 
numerical troubles. Actually this feature can also be appreciated in Figure 3 of Ref. [30|, where the integration starts 
even before recombination. 
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FIG. 25. The deceleration parameter q for the MJW, Starobinsky and Hu-Sawicky models. 
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FIG. 26. The jerk parameter j for the MJW, Starobinsky and Hu-Sawicky models. 



VII. DISCUSSION 

Modified theories of gravity, like metric f(R) theories, have been analyzed thoroughly in recent years mainly to 
explain the accelerated expansion of the Universe which was inferred from the measurements of luminous distance 
of SNIa, while some other theories have been put forward to avoid the need of dark matter. In this article we have 
focused mainly on the mechanism to produce an accelerated expansion. Such behavior can be most easily accounted 
within the general theory of relativity and by (re) introducing the cosmological constant which was proposed by 
Einstein almost one hundred years ago. The motivation behind the proposal for modifying GR was to avoid the 
introduction of such constant and then to circumvent the apparent problems associated with it. Furthermore, such an 
alternative exempts us from adding new fields (scalar or otherwise) like quintessence or k-essence, in order to explain 
that acceleration. Nevertheless this alternative adds, to our opinion much more troubles than solutions. While many 
specific f(R) models are able to produce an accelerated expansion similar to the ACDM model, they have, at the 
same time, spoiled many of the successes of GR or are inconsistent with other features of cosmology, like an adequate 
matter dominated era. Only a few models have succeeded in explaining, at least partially, the actual cosmological 
evolution of the Universe without disturbing, for instance, the predictions at Solar System scales. Since, f(R) theories 
do not introduce a new fundamental principle of nature, there is then, not a deep criteria that favors one among the 
apparently viable f(R) models. Basically some kind of "handcraft" have been used so far to mold specific f(R) or 
in other cases even reconstruction methods [53l |63|. At any rate, simplicity would be in favor of GR with A. In this 
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FIG. 27. Distance modulus for the MJW, Starobinsky and Hu-Sawicky models. The ACDM model is also plotted for reference. 
The data was taken from Riess et al. @]. Here Ho = 70kms~ 1 Mpc~ 1 . 
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FIG. 28. Left: Same as Fig. [27] including the Riess et al. data (green) [j| and the Union 2 compilation of Amanullah et al. 
(red) 0| (color online). Right: similar as the left box but using a different scale and larger z (only Union 2 data is considered 
here) . 




FIG. 29. Similar to Fig. [28] for the luminous distance and larger z (observational data is not included here). 
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FIG. 30. Angular diameter distance for the MJW, Starobinsky and Hu-Sawicky models. The AC'DM model is also plotted for 
reference. 
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FIG. 31. The second derivative f RR := d 2 f/dR 2 of Fig-U as a function of z. Notice that fun goes to zero for large 



article we argued that even if this kind of modified theories can be internally consistent, some special care has to be 
taken into account when they are analyzed. For instance, it was "discovered" that f(R) theories can be equivalent to 
scalar-tensor theories of gravity. This identification is free of inconsistencies provided that the mapping between both 
representations is well defined. For that it is required that the function fa be a monotonic function of the Ricci scalar 
R. Several of the apparently successful models fail to fulfill this condition in general, and yet, different authors have 
used the STT representation. One of the consequences is that the scalar field potential turns to be multivalued. This 
pathology has contributed to create confusion in the subject, in addition to the already existing confusion between 
frames in STT. We have emphasized that in order to avoid such potential drawbacks f(R) theories should be treated 
using the original variables. This is not only possible, but in our view, it turns to be much more transparent, even if 
the equations seem more involved at first sight. Although some other people share this view and have used the original 
variables in several applications, we give a step forward and propose here a general system of equations simpler than 
the one usually used, and when applied particularly to cosmology, it reduces further and can be treated numerically 
as an initial value problem, where the initial values are restricted by a modified Hamiltonian constraint. Previously 
we presented a system of equations that can be used to construct compact objects in static and spherically symmetric 
spacetimes and showed the way to treat them numerically [28]. In the current article we analyzed three specific f(R) 
models that are viable, at least in the background, since they provide an adequate matter dominated behavior followed 
by a correct accelerated expansion. It has been argued that the Miranda et al. model [30| is ultimately incompatible 
when cosmology is analyzed at the level of perturbations or in the Solar System [42[ , but very likely a deeper analysis 
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FIG. 32. Same as Fig 13 II in logarithmic scale. 



is required in order to rule out this model completely [43J. In the future we plan to analyze some other models, like the 
exponential ones [64U66T ] which can be viable as well. For each of the three classes of f(R) models that we considered 
here, we analyzed three possible inequivalent equations of state associated with the modified gravity and which have 
been studied in the past by several authors. Such EOS arise from general EMT that represent the modified gravity or 
geometric dark energy, although such interpretation is to be handled with care as several authors have warned before. 
We have identified at least four inequivalent definitions of the EMT which lead to the three EOS just alluded. Two 
of such EMT's have the unpleasant feature that are not conserved. One of the other four EOS can diverge at some 
red-shift. Thus the more appealing definition comes from Recipe I whose EMT is conserved and the corresponding 
EOS behaves appropriately. We emphasized that it is important to come to an agreement on which definition of the 
EOS will be ultimately compared with the cosmological observations. This is crucial in view that the forthcoming 
data might differentiate between them as precision is increased || Ht| ■ 

It is natural to ask about the new predictions that f(R) theories can make in addition to explaining the accelerated 
expansion while reproducing, in the viable cases, the previous successes of GR and in particular the ACDM model. 
There are indeed several predictions that, in principle, would allow us to distinguish between GR and f(R) theories. 
One of them is the additional polarization modes of gravitational waves like the "breathing" mode [66, 68] . At present 
time the current gravitational- wave detectors are not sensitive enough to detect gravitational radiation and therefore, 
this feature cannot be used as a tool to falsify alternative theories of gravity yet, but in a near future this will be 
certainly possible [69(. The binary pulsar on the other hand can be an excellent tool to do so at present time. We 
know that the binary pulsar can constrain STT even if they succeed in passing the Solar System tests [7(J, thus, 
the same might happen for f(R) theories. At the cosmological lev el, f (R) theories predict a large-scale integrated 
Sachs- Wolf effect which is different from the one predicted in GR [65|, I7M73T ] . This is because the total EMT that 
one can associate with such theories is not necessarily "isotropic" at the time of recombination and thus the metric 
potentials of scalar perturbations in the Newtonian gauge are not equal (in absolute value). Nevertheless, due to the 
so called cosmic variance, it turns difficult to use it as a further constraint for some values of the model parameters. 
On the other hand, weak and strong gravitational lensing as well as the growth of matter perturbations can also be 
affected when those potentials are different [5(J IzE [zH UM, an< ^ so they can further constrain the specific models. 

In this article we have not attempted to best-fitting the parameters using current observations, but rather using 
the KCDM models as a standard, taking a zero spatial curvature as a prior. The former analysis will be pursued in 
a future work. 
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VIII. APPENDIX 



DIMENSIONLESS FORM OF THE COSMOLOGICAL EQUATIONS AND TESTS. 

The differential equations (|50[) ~ (I53I) can be recasted in dimensionless form which is more suitable for numerical 
integration. One can introduce the following dimensionless quantities: 



a = a /do , 


(64) 


t = t/Ho 1 , 


(65) 


H = H/H , 


(66) 


f = f/H 2 , 


(67) 


R — R/H , 


(68) 


Irr — Irr/ H 2 , 


(69) 


Irrr = IrrrJ H 4 , 


(70) 


P = P/Pcrit , 


(71) 


Prad = Prad/Pcrit ' 


(72) 


T = T/p° crit , 


(73) 
(74) 



where p° rit := 3Hq /(8itG). Notice that f R is already dimensionless. 
In terms of these quantities Eqs. (l5P|) - ((54"]) read 



R" 



-r! i 



R \ 1 
6H^J ~ 3f RR H 2 



[z!rrrH 2 R' 2 + 2] - f R R + ?,T] , 



(75) 



H' 



-2H 



R 
6H 



(76) 



H 2 



1 



f RR H 2 R'--(f R R~f) 



fa 



(77) 



(78) 



dt 
da 



1/H , 



(79) 



t — — p and T = — (pbar + Pdm)- Notice now, as compared with Eqs. (|50"j) . ((52^1 . and (|55t . that the 
factor k = 8irG no longer appears, and a factor of '3' appears in Eq. (|75|) next to T in order for T to be given in units 
of Pcj-if Similarly in Eqs. (1771) and (f78|) the factors k/3 have now disappeared for p to be given in units of p° rit . 

As concerns the internal test performed to check the numerical accuracy and consistency of our computer code, we 
checked up to which extent the modified Hamiltonian constraint Eq. (|77l) is verified. We compare the value of H given 
by Eq. ([77| and call it H ana with the value that arises from integrating Eq. (|76p or ((75)1 . that we call it -ff num - The 
relative difference between both measures the degree of verification (or failure) of the constraint. Figure [33] depicts 
this difference as a function of z. Notice that the maximum relative error is lower than 1CP 10 . At the initial z Eq. ((77)) 
is satisfied exactly by construction and thus the precision is infinite. This value is thus not depicted. 

We conclude then that results obtained from the FORTRAN code used to solve the differential equations Eqs. ((75)) 
(1791 and Eq. J62j are reliable. 
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FIG. 33. Relative error commited in the verification of the modified Hamiltonian constraint during the numerical integration 
for the MJW, Starobinsky and Hu-Sawicky models. 
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